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Abstract: Many image processing problems involve identifying the region in the image domain 
occupied by a given entity in the scene. Automatic solution of these problems requires models that 
incorporate significant prior knowledge about the shape of the region. Many methods for including 
such knowledge run into difficulties when the topology of the region is unknown a priori, for example 
when the entity is composed of an unknown number of similar objects. Higher-order active contours 
(HOACs) represent one method for the modelling of non-trivial prior knowledge about shape without 
necessarily constraining region topology, via the inclusion of non-local interactions between region 
boundary points in the energy defining the model. 

The case of an unknown number of circular objects arises in a number of domains, e.g. medical, 
biological, nanotechnological, and remote sensing imagery. Regions composed of an a priori un- 
known number of circles may be referred to as a 'gas of circles'. In this report, we present a HOAC 
model of a 'gas of circles'. In order to guarantee stable circles, we conduct a stability analysis via 
a functional Taylor expansion of the HOAC energy around a circular shape. This analysis fixes one 
of the model parameters in terms of the others and constrains the rest. In conjunction with a suitable 
likelihood energy, we apply the model to the extraction of tree crowns from aerial imagery, and show 
that the new model outperforms other techniques. 

Key-words: tree, crown, extraction, aerial image, higher-order, active contour, gas of circles, prior, 
shape 
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Un modele de contours actifs d'ordre superieur d'un 'gaz de 
cercles' et son application a Fextraction des houppiers 

Resume : Plusieurs problemes en traitement d' image demandent 1' identification de la region dans le 
domaine de 1' image qui correspond a une entite donnee dans la scene. La solution automatique de ces 
problemes demande des modeles qui incorporent une connaissance a priori importante de la forme 
de la region. Beaucoup de methodes pour 1' inclusion d'une telle connaissance ont des difficultes 
quand la topologie est a priori inconnue, par exemple quand 1' entite est composee d'un nombre 
inconnu d'objets similaires. Les contours actifs d'ordre superieur (CAOS) sont une methode pour la 
modelisation d'une connaissance a priori sur la forme non-triviale sans necessairement contraindre 
la topologie de la region, via 1' inclusion d' interactions non-locales entre les points du bord de la 
region dans I'energie qui definie le modele. 

Le cas d'un nombre inconnu d'objets circulaires se souleve dans plusieurs domaines, e.g. I'imagerie 
medicale, biologique, nanotechnologique, et la teledetection. Des regions composees d'un nombre 
inconnu de cercles peuvent etre appelees un "gaz de cercles". Dans ce rapport, nous presentons un 
modele CAOS d'un gaz de cercles. Pour garantir des cercles stables, nous faisons une analyse de sta- 
bilite via un developpement de Taylor fonctionnel autour d'une forme circulaire. Cette analyse fixe 
un des parametres du modele en fonction des autres, et contraint le reste. En combinaison avec une 
energie d' attache aux donnees appropriee, nous appliquons le modele a 1' extraction des houppiers 
dans des images aeriennes, et montrons que sa performance est meilleure que d' autres techniques. 

Mots-cles : arbre, houppier, extraction, image aerienne, ordre superieur, contour actif, gaz de 
cercles, a priori, forme 
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1 Introduction 

Forestry is a domain in which image processing and computer vision techniques can have a signifi- 
cant impact. Resource management and conservation, whether commercial or in the pubHc domain, 
require information about the current state of a forest or plantation. Much of this information can 
be summarized in statistics related to the size and placement of individual tree crowns {e.g. mean 
crown area and diameter, density of the trees). Currently, this information is gathered using ex- 
pensive field surveys and time-consuming semi-automatic procedures, with the result that partial 
information from a number of chosen sites frequently has to be extrapolated. An image processing 
method capable of automatically extracting tree crowns from high resolution aerial or satellite im- 
ages (an example is shown in figure 1) and computing statistics based on the results would greatly 
aid this domain. 



The tree crown extraction problem can be viewed as a special case of a general image under- 
standing problem: the identification of the region R in the image domain corresponding to some 
entity or entities in the scene. In order to solve this problem in any particular case, we have to 
construct, even if only implicitly, a probability distribution on the space of regions P(i^|/, K). This 
distribution depends on the current image data / and on any prior knowledge K we may have about 
the region or about its relation to the image data, as encoded in the likelihood P(/|i^, K) and the 
prior P(i? I iiT) appearing in the Bayes' decomposition of P(i?|/, K) (or equivalently in their energies 




Figure 1 : Real image with planted forest ©IFN. 
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— hiV{I\R^K) and — hiV{R\K)). This probability distribution can then be used to make estimates 
of the region we are looking for. 

In the automatic solution of realistic problems, the prior knowledge K, and in particular prior 
knowledge about the 'shape' of the region, as described by V{R\K), is critical. The tree crown 
extraction problem provides a good example: particularly in plantations, R takes the form of a 
collection of approximately circular connected components of similar size. There is thus a great 
deal of prior knowledge about the region sought. The question is then how to incorporate such prior 
knowledge into a model for Rl If the model does not include enough prior knowledge, it will be 
necessary to include it in some other form. So, for example, the use of classical active contour 
energies to find entities in images usually requires the initialization of the contour close to the entity 
to be found, which represents a large injection of prior knowledge by the user. 

The simplest prior information concerns the smoothness properties of the boundary of the re- 
gion. For example, the Ising model and many active contour models [3, 4, 7, 8, 18] use the length 
of the region boundary and its interior area as their prior energies. This type of prior information 
can be augmented using more demanding measures of smoothness, for example the boundary curva- 
ture [14, 18]. These models are all limited, though, by the fact that they are integrals over the region 
boundary of some function of various derivatives of the boundary. In consequence, they capture lo- 
cal differential geometric information, corresponding to local interactions between boundary points, 
but can say nothing more global about the shape of the region. 

To go further, it is therefore clear that one must introduce longer range interactions. There 
are two principal ways to do this: one is to introduce hidden variables, given which the original 
variables of interest are (more or less) independent. Marginalizing over the hidden variables then 
introduces interactions between the original variables. Another is to include explicit long-range 
interactions between the original variables (these interactions may also have an interpretation in 
terms of marginalization over some hidden variables, which are left implicit). 

The first approach has been much investigated, in the form of template shapes and their defor- 
mations [5, 9, 10, 11, 12, 13, 17, 21, 22, 23, 24, 26]. Here a probability distribution or an energy 
is defined based on a distance measure of some kind between regions. One region, the template, is 
fixed, while the other is the variable R. This type of model constrains R to be close to the template 
region in the space of regions. Template regions may be learned from examples or fixed by hand; 
similarly the distance function maybe based, for example, on the learned covariance of a Gaussian 
distribution, or chosen a priori. The most sophisticated methods use the kernel trick to define the 
distance as a pullback from a high-dimensional space, thereby allowing more complex behaviours. 
Multiple templates may also be used, corresponding to a mixture model. It is clear that these meth- 
ods implicitly introduce long-range interactions: if you know that one half of a given region aligns 
well with the template, this tells you something about the likely position and shape of the other half. 

The above methods assign high probability to regions that lie 'close' to certain points in the 
space of regions. As such, it is difficult to construct models of this type that favour regions for which 
the topology, and in particular the number of connected components, is unknown a priori. There 
are many problems, however, for which this is the case, for example, the extraction of networks, or 
the extraction of an unknown number of objects of a particular type from astronomical, biological. 
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medical, or remote sensing images. For this type of prior knowledge, a different type of model is 
needed. Higher-order active contours (HOACs) are one such category of models. 

HOACs, first described by Rochery et al. [ >] (see also [30, 31] for fuller descriptions), take 
the second approach mentioned above. They introduce explicit long-range interactions between 
region boundary points via energies that contain multiple integrals over the boundary, thus avoiding 
the use of template shapes. HOAC energies can be made intrinsically Euclidean invariant, and, 
as required by the above analysis, incorporate sophisticated prior information about region shape 
without necessarily constraining the topology of the region. As with other methods incorporating 
significant prior knowledge, it is not necessary to introduce extra knowledge via an initialization 
close to the target region: a generic initialization suffices, thus rendering the method quasi- automatic. 
Rochery et al. [ ] applied the method to road extraction from satellite and aerial images using a prior 
which favours network- like objects. 

In this report, we describe a HOAC model of a 'gas of circles': the model favours regions 
composed of an a priori unknown number of circles of a certain radius. For such a model to work, 
the circles must be stable to small perturbations of their boundaries, i.e. they must be local minima 
of the HOAC energy, for otherwise a circle would tend to 'decay' into other shapes. This is a non- 
trivial requirement. We impose it by performing a functional Taylor expansion of the HOAC energy 
around a circle, and then demanding that the first order term be zero for all perturbations, and that 
the second order term be positive semi-definite. These conditions allow us to fix one of the model 
parameters in terms of the others, and constrain the rest. Experiments using the HOAC energy 
demonstrate empirically the coherence between these theoretical considerations and the gradient 
descent algorithm used in practice to minimize the energy. 

The model has many potential applications, to medical, biological, physical, and remote sensing 
imagery in which the entities to be identified are circular. We choose to apply it to the tree crown 
extraction problem from aerial imagery, using the 'gas of circles' model as a prior energy, and an 
appropriate likelihood. We will see that the extra prior knowledge included in the 'gas of circles' 
model permits the separation of trees that cannot be separated by simpler methods, such as maximum 
likelihood or classical active contours. 

In the rest of this section, we present a brief introduction to HOACs. In section 2, we describe 
the 'gas of circles' HOAC model. The key to this model is the analysis of the stability of a circle 
as a function of the model parameters. To demonstrate the prior knowledge contained in the model, 
and the empirical correctness of the stability analysis, we present experimental results using the new 
energy. In section 3, we apply the new model to tree crown extraction. We describe a likelihood 
energy for trees based on the image intensity and gradient, and then present experimental results on 
synthetic data and on aerial images. We conclude in section 4, and discuss some open issues with 
the model. 

1.1 Higher order active contours 

As described in section 1 , HOACs introduce long-range interactions between boundary points not 
via the intermediary of a template region or regions to which R is compared, but directly, by using 
energy terms that involve multiple integrals over the boundary. The integrands of such integrals thus 
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depend on two or more, perhaps widely separated, boundary points simultaneously, and can thereby 
impose relations between tuples of points. Euclidean invariance of such energies can be imposed 
directly on the energy, without the necessity to estimate a transformation between the boundary 
sought and the template. Perhaps more importantly, because there is no template, the topology of 
the region need not be constrained, a factor that is critical when the topology is not known a priori. 

As with all active contour models, a region R is represented by its boundary, dR. There are 
various ways to think of the boundary of a region. If the region has only one connected compo- 
nent, which is also simply-connected, then a boundary is an equivalence class of embeddings of the 
circle under the action of orientation-preserving diffeomorphisms of . When more, possibly 
multiply-connected components are included, however, things get complicated. First, the number of 
embeddings of S*^ that are required depends on the topology, and second, there are constraints on 
the orientations of different components if they are to represent regions with handles. 

An alternative is to view dR as a closed 1 -chain 7 in the image domain 1] ([ ] is a useful 
reference for the following discussion). Although region boundaries correspond to a special subset 
of closed 1 -chains known as domains of integration, active contour energies themselves are defined 
for general 1 -chains. It is convenient to use this more general context to distinguish HOAC energies 
from classical active contours, because it allows for notions of linearity to be used to characterize 
the complexity of energy functionals. 

Using this representation, HOAC energies can be defined as follows [30, 31]. Let 7 be a 1-chain 
in Vt, and ^7 be its domain. Then 7^ : (□7)^ is an n-chain in QJ^ . We define a class of 

(n — p)-forms on Vf^ that are 1-forms with respect to (n — p) factors and 0-forms with respect to the 
remaining p factors (by symmetry, it does not matter which p factors). These forms can be pulled 
back to (□7)^ by 7"^. The Hodge duals of the p 0-form factors with respect to the induced metric 
on ^7 can then be taken independently on each such factor, thus converting them to 1-forms, and 
rendering the whole form an n-form on (□7)^. This n-form can then be integrated on (□7)^. 

In the (n,p) = (n, 0) cases, we are simply integrating a general n-form on the image of 7^ 
in 1]^, thus defining a linear functional on the space of n-chains in OP', and hence an n*-order 
monomial on the space of 1 -chains in ^. Taking arbitrary linear combinations of such monomials 
then gives the space of polynomial functionals on the space of 1 -chains. By analogy we refer to 
the general (n,p) cases as 'generalized n*-order monomials' on the space of 1-chains in Vt, and 
to arbitrary linear combinations of the latter as 'generalized polynomial functionals' on the space 
of 1-chains in Vt. HOAC energies are generalized polynomial functionals. Standard active contour 
energies are generalized linear functionals on 1-chains in this sense, hence the term 'higher-order'. 

The (1,1) case is simply the boundary length in some metric. An interesting application of the 
(2,2) case to topology preservation is described by Sundaramoorthi and Yezzi [ ]. The (1,0) case 
gives the region area in some metric. 

To be more concrete, we specialize to the (n, 0) case. Let F be an n-form on . We pull F 
back to the domain of 7"^ and integrate it: 




(1.1) 



RR n° 0123456789 



8 



Horvdth & Jermyn & Kato & Zerubia 



Specializing again to the case n = 2, and using the antisymmetry of F together with the symmetry 
of 7^, we can rewrite the energy functional in this case as 



^(7) = / F= (7 X 7)*F = / / dt dt' T{t) • F(7(t), 7(^0) • r{t') , (1.2) 

where F{x^x'), for each {x^x') G is a 2 x 2 matrix, t is a coordinate on ^7, and r = 7 is the 
tangent vector to 7. 

By imposing Euclidean invariance on this term, and adding linear terms, Rochery et al. [ ] 
defined the following higher-order active contour prior: 

Eg(7) = \cLin) + acAin) jj dtdt' Tit') • r(t) ^{R{t, t')) , (1.3) 

where L is the boundary length functional, A is the interior area functional and R{t^ t') = \j{t) — 
7(t')| is the Euclidean distance between j{t) and 7(t'). Rochery et al. [ ] used the following 
interaction function <l>: 

z < d — e ^ 

. £^ _ 1 sin ZL(£^) d-e<z <d^e , (1.4) 
z>d^e. 

In this paper, we use this same interaction function with d = e, but other monotonically decreas- 
ing functions lead to qualitatively similar results. 




2 The ^gas of circles' model 

For certain ranges of the parameters involved, the energy (1.3) favours regions in the form of net- 
works, consisting of long narrow arms with approximately parallel sides, joined together at junctions, 
as described by Rochery et al. [29, 30, 31]. It thus provides a good prior for network extraction from 
images. This behaviour does not persist for all parameter values, however, and we will exploit 
this parameter dependence to create a model for a 'gas of circles', an energy that favours regions 
composed of an a priori unknown number of circles of a certain radius. 

For this to work, a circle of the given radius, hereafter denoted ro, must be stable, that is, it 
must be a local minimum of the energy. In section 2.1, we conduct a stability analysis of a circle, 
and discover that stable circles are indeed possible provided certain constraints are placed on the 
parameters. More specifically, we expand the energy Eg in a functional Taylor series to second order 
around a circle of radius ro . The constraint that the circle be an energy extremum then requires that 
the first order term be zero, while the constraint that it be a minimum requires that the operator in the 
second order term be positive semi-definite. These requirements constrain the parameter values. In 
subsection 2.2, we present numerical experiments using Eg that confirm the results of this analysis. 
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2.1 Stability analysis 

We want to expand the energy Eg around a circle of radius tq . We denote a member of the equiva- 
lence class of maps representing the 1 -chain defining the circle by 70. The energy Eg is invariant to 
diffeomorphisms of 070, and thus is well-defined on 1 -chains. To second order, 

5E 1 6"^ E 

E,ij) = E.i^o + <57) = E.ijo) + ihl^ho + ^(^^I^IS^K ■ (2-1) 

where (•!•) is a metric on the space of 1 -chains. 

Since 70 represents a circle, it is easiest to express it in terms of polar coordinates r, 6> on Vt. 
For a suitable choice of coordinate on 5^, a circle of radius ro centred on the origin is then given 
by 7o(^) = (^o(^), ^o(^)), where ro(t) = ro, 0{t) = t, and t G [— 7r,7r). We are interested in the 
behaviour of small perturbations (^7 = ((5r, 59). The first thing to notice is that because the energy 
Eg is defined on 1 -chains, tangential changes in 7 do not affect its value. We can therefore set 
(56> = 0, and concentrate on 8r. 

On the circle, using the arc length parameterization t, the integrands of the different terms in 
Eg are functions oft — t' only; they are invariant to translations around the circle. In consequence, 
the second derivative S'^Eg/Sj{t)Sj{t') is also translation invariant, and this implies that it can be 
diagonalized in the Fourier basis of the tangent space at 70. It thus turns out to be easiest to perform 
the calculation by expressing Sr in terms of this basis: 

Sr{t) = ^ake^-^'\ (2.2) 

k 

where k G {m/ro : m G Z}. Below, we simply state the resulting expansions to second order in 
the a/c for the three terms appearing in equation (1.3). Details can be found in appendix A. 
The boundary length and interior area of the region are given to second order by 




(2.3) 



A(7)= / dO dr' r' = iTrl^2iTrQaQ^iTYjc^k? ' (2-4) 

Note the /c^ in the second order term for L. This is the same frequency dependence as the Laplacian, 
and shows that the length term plays a similar smoothing role for boundary perturbations as the 
Laplacian does for functions. In the area term, by contrast, the Fourier perturbations are 'white 
noise'. 

It is also worth noting that there are no stable solutions using these terms alone. For the circle to 
be an extremum, we require Ac27r + ac^nvQ = 0, which tells us that ac = — Ac/^o- The criterion 
for a minimum is, for each k, Xcr^k'^ + ^ 0. Note that we must have Ac > for stability at 
high frequencies. Substituting for ac, the condition becomes Xc{rok'^ ~ ^0^) ^ 0- Substituting 
k = m/ro, gives the condition — 1 > 0. Two points are worth noting. The first is the one we 
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have already made: the zero frequency perturbation is not stable. The second is that the m = 1 
perturbation is marginally stable to second order, that is, such changes require no energy to this 
order. To fully analyse them, we must therefore go to higher order in the Taylor series. This feature 
will appear also in the analysis of the full energy E^. 

The quadratic term can be expressed to second order as 

n p7T pit ptt 

/ / dt dt' T{t') ■ T{t) t')) = 2tt dp Foo(p) + inao / dp Fio{p) 

J J —TV J —TT J —T^ 

+ ^27r|afe|2|[2 f dp F2o{p) + T dp 6"''^°'^^ F21 (p)" 

- \liT^k dpe-^-^^^F2s{p)] + [rlk^ f dpe-^^^^^FM] | , (2.5) 

The Fij are functionals of ^ (hence functions of d and e for <l> given by equation (1 .4)), and functions 
of ro, as well as functions of the dummy variable p. 

Combining equations (2.3), (2.4), and (2.5), we find the energy functional (1.3) up to the second 
order: 

Egijo^^l) = eo(ro) +aoei(ro) + ^^\ak\^e2{k,ro) 

^ k 

= |27rAcro + Tracr^ - 7rf3cGoo{ro)^ + ao|27rAc + ^nacro - 27r/3cGio(ro)| 
+ I ^\akf{'27rXcrok^ + 27rac 

k 

- 27rPc [2G2o(ro) + G2i{k, ro) - 2irokG23{k, ro) + rieG24{k, ro)] } , 



(2.6) 

where Gij = J^^ dp e~'^^^^^~^^^^^^^ Fij (p). Note that as anticipated, there are no off-diagonal terms 
linking ak and a^/ for k ^ k': the Fourier basis diagonalizes the second order term. 

2.1.1 Parameter constraints 

Note that a circle of any radius is always an extremum for non-zero frequency perturbations (ak for 
/c 7^ 0), as these Fourier coefficients do not appear in the first order term (this is also a consequence 
of invariance to translations around the circle). The condition that a circle be an extremum for ao as 
well (ei = 0) gives rise to a relation between the parameters: 

Pc[Xc, ac, ro) = ^ . , (2.7) 

where we have introduced ro to indicate the radius at which there is an extremum, to distinguish 
it from ro, the radius of the circle about which we are calculating the expansion (2.1). The left 
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hand side of figure 2 shows a typical plot of the energy eo of a circle versus its radius ro, with 
the Pc parameter fixed using the equation (2.7) with Ac = 1.0, a = 0.8, and fo = 1-0. The 
energy has a minimum at ro = ro as desired. The relationship between fo and Pc is not quite as 
straightforward as it might seem though. As can be seen, the energy also has a maximum at some 
radius. It is not a priori clear whether it will be the maximum or the minimum that appears at fo. 
If we graph the positions of the extrema of the energy of a circle against (3c for fixed ac, we find a 
curve qualitatively similar to that shown in figure 3 (this is an example of a fold catastrophe). The 
solid curve represents the minimum, the dashed the maximum. Note that there is indeed a unique 
f3c for a given choice of f o . Denote the point at the bottom of the curve by (p'^c^^o^)' Note that 
at Pc = Pc^ , the extrema merge and for Pc < Pq^ , there are no extrema: the energy curve is 
monotonic because the quadratic term is not strong enough to overcome the shrinking effect of the 
length and area terms. Note also that the minimum cannot move below ro = r^^ . This behaviour is 
easily understood qualitatively in terms of the interaction function in equation (1.4). If 2ro < d — e, 
the quadratic term will be constant, and no force will exist to stabilize the circle. In order to use 
equation (2.7) then, we have to ensure that we are on the upper branch of figure 3. 




0.5 1 1.5 2 2.5 3 1 2 3 4, 5 6 7 8 

Figure 2: Plots of eo against ro and 62 against rok. Left: the energy of a circle eo plotted against 
radius ro for Ac = 1.0, a = 0.8, and Pc = 1.39 calculated from equation (2.7) with fo = 1.0. (The 
parameters of $ are (i = 1.0 and e = 1.0, but note that it is not necessary in general that d = ro.) 
The function has a minimum at ro = fo as desired. Right: the second derivative of Eg, 62, plotted 
against r^k for the same parameter values. The function is non-negative for all frequencies. 

Equation (2.7) gives the value of Pc that provides an extremum of eo with respect to changes 
of radius ao at a given fo (ei(fo) = 0), but we still need to check that the circle of radius fo is 
indeed stable to perturbations with non-zero frequency, i.e. that e2(/c, fo) is non-negative for all k. 
Scaling arguments mean that in fact the sign of e2 depends only on the combinations = r^/d and 
Oic = {d/\c)oLc- The equation for e2 can then be used to obtain bounds on ac in terms of fo. 
(Details of these calculations and bounds will be given elsewhere.) The right hand side of figure 2 



RR n° 0123456789 



12 



Horvdth & Jermyn & Kato & Zerubia 



shows a plot of 62 (/c, fo) against rok for the same parameter values used for the right hand side, 
showing that it is non-negative for all rok. 




Figure 3: Schematic plot of the positions of the extrema of the energy of a circle versus (3c- 

We call the resulting model, the energy Eg with parameters chosen according to the above crite- 
ria, the 'gas of circles' model. 

2.2 Geometric experiments 

In order to illustrate the behaviour of the prior energy with parameter values fixed according to 
the above analysis, in this section we show the results of some experiments using this energy (there 
are no image terms). Figure 4 shows the result of gradient descent using Eg starting from various 
different initial regions. (For details of the implementation of gradient descent for higher-order 
active contour energies using level set methods, see [30, 3 1 ] .) In the first column, four different initial 
regions are shown. The other three columns show the final regions, at convergence, for three different 
sets of parameters. In particular, the three columns have fo = 15.0, 10.0, and 5.0 respectively. 

In the first row, the initial shape is a circle of radius 32 pixels. The stable states, which can be 
seen in the other three columns, are circles with the desired radii in every case. In the second row, 
the initial region is composed of four circles of different radii. Depending on the value of fo, some 
of these circles shrink and disappear. This behaviour can be explained by looking at figure 2. As 
already noted, the energy of a circle eo has a maximum at some radius rmax- If an initial circle has a 
radius less than rmax, it will 'slide down the energy slope' towards ro = 0, and disappear. If its radius 
is larger than rmax, it will finish in the minimum, with radius fo. This is precisely what is observed 
in this second experiment. In the third row, the initial condition is composed of four squares. The 
squares evolve to circles of the appropriate radii. The fourth row has an initial condition composed 
of four differing shapes. The nature of the stable states depends on the relation between the stable 
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radius, fo, and the size of the initial shapes. If fo is much smaller than an initial shape, this shape 
will 'decay' into several circles of radius fo. 




(Initial) (fo = 15) (fo = 10) (fo = 5) 



Figure 4: Experimental results using the geometric term: the first column shows the initial condi- 
tions; the other columns show the stable states for various choices of the radius. 



3 Data terms and experiments 

In this section, we apply the 'gas of circles' model developed in section 2 to the extraction of trees 
from aerial images. This is just one of many possible applications, corresponding to the mission of 
the Ariana research group. In the next section, we give a brief state of the art for tree crown extrac- 
tion, and then present the data terms we use in section 3.2. In section 3.3, we describe tree crown 
extraction experiments on aerial images and compare the results to those found using a classical ac- 
tive contour model. In section 3.4, we examine the robustness of the method to noise using synthetic 
images. This illuminates the principal failure modes of the model, which will be further discussed 
in section 4, and which point the way for future work. In section 3.5, we illustrate the importance 
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of prior information via tree crown separation experiments on synthetic images, and compare the 
results to those obtained using a classical active contour model. 

3.1 Previous work 

The problem of locating, counting, or delineating individual trees in high resolution aerial images 
has been studied in a number of papers. 

Gougeon [ 16] observes that trees are brighter than the areas separating them. Local minima of 
the image are found using a 3 x 3 filter, and the 'valleys' connecting them are then found using a 5 x 5 
filter. The tree crowns are subsequently delineated using a five-level rule-based method designed to 
find circular shapes, but with some small variations permitted. Larsen [ , ] concentrates on 
spruce tree detection using a template matching method. The main difference between these two 
papers is the use of multiple templates in the second. The 3D shape of the tree is modelled using 
a generalized ellipsoid, while illumination is modelled using the position of the sun and a clear- 
sky model. Reflectance is modelled using single reflections, with the branches and needles acting 
as scatterers, while the ground is treated as a Lambertian surface. Template matching is used to 
calculate a correlation measure between the tree image predicted by the model and the image data. 
The local maxima of this measure are treated as tree candidates, and various strategies are then used 
to eliminate false positives. Brandtberg and Walter [ ] decompose an image into multiple scales, and 
then define tree crown boundary candidates at each scale as zero crossings with convex greyscale 
curvature. Edge segment centres of curvature are then used to construct a candidate tree crown 
region at each scale. These are then combined over different scales and a final tree crown region is 
grown. Andersen et al. [ ] use a morphological approach combined with a top-hat transformation 
for the segmentation of individual trees. 

All of these methods use multiple steps rather than a unified model. Closer in spirit to the present 
work is that of Perrin et al. [27, 28], who model the collection of tree crowns by a marked point 
process, where the marks are circles or ellipses. An energy is defined that penalizes, for example, 
overlapping shapes, and controls the parameters of the individual shapes. Reversible Jump MCMC 
and simulated annealing are used to estimate the tree crown configuration. Compared to the work 
described in this paper, the method has the advantage that overlapping trees can be represented as 
two separate objects, but the disadvantage that the tree crowns are not precisely delineated due to 
the small number of degrees of freedom for each mark. 

3.2 Data terms and gradient descent 

In order to couple the region model to image data, we need a likelihood, V{I\R^K). The images 
we use for the experiments are coloured infrared (CIR) images. Originally they are composed of 
three bands, corresponding roughly to green, red, and near infrared (NIR). Analysis of the one- 
point statistics of the image in the region corresponding to trees and the image in the background, 
shows that the 'colour' information does not add a great deal of discriminating power compared to a 
'greyscale' combination of the three bands, or indeed the NIR band on its own. We therefore model 
the latter. 
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The images have a resolution ~ O.Sm/pixel, and tree crowns have diameters of the order of ten 
pixels. Very Httle if any dependence remains between the pixels at this resolution, which means, 
when combined with the paucity of statistics within each tree crown, that pixel dependencies {i.e. 
texture) are very hard to use for modelling purposes. We therefore model the interior of tree crowns 
using a Gaussian distribution with mean ji and covariance cf'^Sr, where 5a is the identity operator 
on images on A c Vt. 

The background is very varied, and thus hard to model in a precise way. We use a Gaussian 
distribution with mean Jl and variance In general, ji > Jl, and a < a; trees are brighter and 

more constant in intensity than the background. The boundary of each tree crown has significant 
inward-pointing image gradient, and although the Gaussian models should in principle take care of 
this, we have found in practice that it is useful to add a gradient term to the likelihood energy. Our 
likelihood thus has three factors: 

P(/|i?, K) = gnilR) gjiiiji) foRihn) . 

where Ir and Ir are the images restricted to R and R respectively, and qr and Qr are proportional 
to the Gaussian distributions already described, i.e. 

- \ngR{lR) = d^x ^{Ir{^) - m)' (3.1) 

and similarly for g^. The function Jqr depends on the gradient of the image dl on the boundary 
dR: 

- \nfoR{lQR) = \ f dt n{t) . dl{t) (3.2) 

where n is the unnormalized outward normal to 7. The normalization constant Z is thus a function 
of /i, cr, Jl, a, and A^. Z is also a functional of the region R. To a first approximation, it is a linear 
combination of L{dR) and A{R). It thus has the effect of changing the parameters Ac and ac in Eg. 
However, since these parameters are essentially fixed by hand (the criteria described in section 2.1.1 
only allow us to fix (3c and constrain ac), knowledge of the normalization constant does not change 
their values, and we ignore it once the likelihood parameters have been learnt. 
The full model is then given by E{R) = Ei{I, R) + Eg{R), where 

E,{I,R) = - \nF{I\R,K) ^\nZ = - \ngR{lR) - \ng^{I^) - ln/ai^(/ai?) • 

The energy is minimized by gradient descent. The gradient descent equation for E is 

, (/(7W)-/i)' 



ds^ ' ''''' 2cr2 2^2 



- Xci^it) -ac^Pc I dt' R{t, t') ■ n{t')^{R{t, t')) , (3.3) 
Jn7 

where s is the descent time parameter, ^(t, t') = (7(t) — "){t')) /\^{t) — \ and n is the signed 
boundary curvature. As already mentioned, to evolve the region we use the level set framework of 
Osher and Sethian [25] extended to the demands of nonlocal forces such as equation (3.3) [30, 31]. 
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3.3 Tree crown extraction from aerial images 

In this section, we present the results of the appHcation of the above model to 50 cm/pixel colour 
infrared aerial images of poplar stands located in the 'Saone et Loire' region in France. The images 
were provided by the French National Forest Inventory (IFN). As stated in section 3.2, we model 
only the NIR band of these images, as adding the other two bands does not increase discriminating 
power. The tree crowns in the images are ~ 8-10 pixels in diameter, i.e. ^ 4-5m. 

In the experiments, we compare our model to a classical active contour model (f3c = 0). The 
parameters /i, a, p, and a were the same for both models, and were learned from hand-labelled 
examples in advance. The classical active contour prior model thus has three free parameters (A^, 
Ac and ac), while the HOAC 'gas of circles' model has six (A^, Xc, (^c^ Pc^ d and ro). We fixed 
ro based on our prior knowledge of tree crown size in the images, and d was then set equal to tq. 
Once ac and Ac have been fixed, (3c is determined by equation (2.7). There are thus three effective 
parameters for the HOAC model. In the absence of any method to learn A^, ac and Ac, they were 
fixed by hand to give the best results, as with most applications of active contour models. The values 
of Ai, Ac and ac were not the same for the classical active contour and HOAC models; they were 
chosen to give the best possible result for each model separately. 

The initial region in all real experiments was a rounded rectangle slightly bigger than the image 
domain. The image values in the region exterior to the image domain were set to fl to ensure that the 
region would shrink inwards. 

Figure 5 illustrates the first experiment. On the left is the data, showing a regularly planted poplar 
stand. The result is shown on the right. We have applied the algorithm only in the central part of the 
image, for reasons that will be explained in section 4. 

Figure 6 illustrates a second experiment. On the left is the data. The image shows a small piece 
of an irregularly planted poplar forest. The image is difficult because the intensities of the crowns 
are varied and the gradients are blurred. In the middle is the best result we could obtain using a 
classical active contour. On the right is the result we obtain with the HOAC 'gas of circles' model. ^ 
Note that in the classical active contour result several trees that are in reality separate are merged 
into single connected components, and the shapes of trees are often rather distorted, whereas the 
prior geometric knowledge included when P ^ allows the separation of almost all the trees and 
the regularization of their shapes. 

Figure 7 illustrates a third experiment. Again the data is on the left, the best result obtained with 
a classical active contour model is in the middle, and the result with the HOAC 'gas of circles' model 
is on the right. The trees are closer together than in the previous experiment. Using the classical 
active contour, the result is that the tree crown boundaries touch in the majority of cases, despite 
their separation in the image. Many of the connected components are malformed due to background 
features. The HOAC model produces more clearly delineated tree crowns, but there are still some 
joined trees. We will discuss this further in section 4 

^Unless otherwise specified, in the figure captions the values of the parameters learned from the image are shown when the 
data is mentioned, in the form (fi, a, p., a). The other parameter values are shown when each result is mentioned, in the form 
(Xi, Ac, ac, (^c, d,ro), truncated if the parameters are not present. All parameter values are truncated to two significant 
figures. Unless otherwise specified, images were scaled to take values in [0, 1]. The region boundary is shown in white. 
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Figure 8 shows a fourth experiment. The data is on the left, the best result obtained with a 
classical active contour model is in the middle, and the result with the HOAC 'gas of circles' model 
is on the right. Again, the 'gas of circles' model better delineates the tree crowns and separates more 
trees, but some joined trees remain also. The HOAC model selects only objects of the size chosen, 
so that false positives involving small objects do not occur. 




Figure 5: Left: real image with a planted forest ©IFN (0.3, 0.06, 0.05, 0.05). Right: the result 
obtained using the 'gas of circles' model (529, 5.88, 5.88, 5.64, 4, 4). 




Figure 6: From left to right: image of poplars ©IFN (0.73, 0.11, 0.23, 0.094); the best 
result with a classical active contour (880,13,73); result with the 'gas of circles' model 
(100,6.7,39,31,4.2,4.2). 

Table 1 shows the percentages of correct tree detections, false positives and false negatives (two 
joined trees count as one false negative), obtained with the classical active contour model and the 
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Figure 7: From left to right: image of poplars ©IFN (0.71, 0.075, 0.18, 0.075); the best re- 
sult with a classical active contour (24000,100,500); result with the 'gas of circles' model 
(1500,25,130,100,3.5,3.5). 




Figure 8: From left to right: image of poplars ©IFN (0.71, 0.075, 0.18, 0.075); the best re- 
sult with a classical active contour (35000,100,500); result with the 'gas of circles' model 
(1200,20,100,82,3.5,3.5). 



'gas of circles' model in the experiments shown in figures 6, 7, and 8. The 'gas of circles' model 
outperforms the classical active contour in all measures, except in the number of false negatives in 
the experiment in figure 7. 

Once the segmentation result has been obtained, it is a relatively simple matter to compute statis- 
tics of interest to the forestry industry: number of trees, total area, number and area density, and so 
on. 
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Model 




CAC 






HOAC 




Figure 


CD% 


FP% 


FN% 


CD% 


F+% 


F- % 


Figure 6 


85 





15 


97 





3 


Figure 7 


96.2 


2.8 


1.9 


97.7 





2.3 


Figure 8 


89.4 


5 


5.6 


95.5 


0.6 


3.9 



Table 1 : Results on real images using a classical active contour model (CAC) and the 'gas of circles' 
model (HOAC). CD: correct detections; FP: false positives; FN: false negatives (two joined trees 
count as one false negative). 

3.4 Noisy synthetic images 

In this section, we present the results of tests of the sensitivity of the model to noise in the image. 
Fifty synthetic images were created, each with ten circles with radius 8 pixels and ten circles with 
radius 3.5 pixels, placed at random but with overlaps rejected. Six different levels of white Gaussian 
noise, with image variance to noise power ratios from —5 dB to 20 dB, were then added to the 
images to generate 300 noisy images. Six of these, corresponding to noisy versions of the same 
original image, were used to learn /i, a, /i, and a. The model used was the same as that used for the 
aerial images, except that was set equal to zero. The parameters were adjusted to give a stable 
radius of 8 pixels. 

The results obtained on the noisy versions of one of the fifty images are shown in figure 9. 
Table 2 shows the proportion of false negative and false positive circle detections with respect to the 
total number of potentially correctly detectable circles (500 = 50 x 10), as well as the proportion of 
'joined circles', when two circles are grouped together (an example can be seen in the bottom right 
image of figure 9). Detections of one of the smaller circles (which only occurred a few times even 
at the highest noise level) were counted as false positives. The method is very robust with respect to 
all but the highest levels of noise. The first errors occur at 5 dB, where there is a 2% false positive 
rate. At dB, the error rate is ~ 10%, i.e. one of the ten circles in each image was misidentified on 
average. At —5 dB, the total error rate increases to ~ 30%, rendering the method not very useful. 

Note that the principal error modes of the model are false positives and joined circles. There are 
good reasons why these two types of error dominate. We will discuss them further in section 4. 

3.5 Circle separation: comparison to classical active contours 

In a final experiment, we simulated one of the most important causes of error in tree crown extraction, 
and examined the response of classical active contour and HOAC models to this situation. The errors, 
which involve joined circles similar to those found in the previous experiment, are caused by the fact 
that in many cases nearby tree crowns in an image are connected by regions of significant intensity 
with significant gradient with respect to the background, thus forming a dumbbell shape. Calling 
the bulbous extremities, the 'bells', and the join between them, the 'bar', the situation arises when 
the bells are brighter than the bar, while the bar is in turn brighter than the background, and most 
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* * * 




(0.90, 0.028, 0.11, 0.028, 100, 100, 170) (0.85, 0.043, 0.16, 0.043, 33, 33, 58) 





(0.78, 0.061, 0.23, 0.062, 13, 13, 22) (0.71, 0.081, 0.31, 0.081, 4, 4, 6.9) 





(0.65, 0.098, 0.37, 0.099, 1.4, 1.4, 2.5) (0.60, 0.11, 0.43, 0.11, 0.51, 0.51, 0.89) 



Figure 9: One of the synthesized images, with six different levels of added white Gaussian noise. 
Reading from left to right, top to bottom, the image variance to noise power ratios are 20, 15, 10, 
5, 0, —5 dB. Parameter values in the form (/i, a, /i, a, Ac, cec, f^c) are shown under the six images. 
The parameters d and tq were fixed to 8 throughout. 
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15 











10 
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6.4 


4 





-5 


27.6 


3.6 


23 



Table 2: Results on synthetic noisy images. FP, FN, J: percentages of false positive, false nega- 
tive, and joined circle detections respectively, with respect to the potential total number of correct 
detections. 

importantly, the gradient between the background and the bar is greater than that between the bar 
and the bells. 

The first row of figure 10 shows a sequence of bells connected by bars. The intensity of the 
bar varies along the sequence, resulting in different gradient values. We applied the classical active 
contour and HOAC 'gas of circles' models to these images. 

The middle row of figure 10 shows the best results obtained using the classical active contour 
model. The model was either unable to separate the individual circles, or the region completely 
vanished. The intuition is that if there is insufficient gradient to stop the region at the sides of the 
bar, then there will also be insufficient gradient to stop the region at the boundary between the bar 
and the bells, so that the region will vanish. On the other hand, if there is sufficient gradient between 
the bar and the background to stop the region, the circles will not be separated, and a 'bridge' will 
remain between the two circles.^ 

The corresponding results using the HOAC 'gas of circles' model can be seen in the bottom row 
of figure 10. All the circles were segmented correctly, independent of the gray level of the junction. 
Encouraging as this is, it is not the whole story, as we indicated in section 3.4. We make a further 
conmient on this issue in section 4, which now follows. 

4 Conclusion 

Higher-order active contours allow the inclusion of sophisticated prior information in active contour 
models. This information can concern the relation between a region and the data, i.e. the likelihood 
V{I\R^K), but more often it concerns the prior probability V{R\K) of a region, or in other words, 
its 'shape'. HOACs are particularly well adapted to including shape information about regions for 
which the topology is unknown a priori. 

In this paper, we have shown via a stability analysis that a HOAC energy can be constructed that 
describes a 'gas of circles', that is, it favours regions composed of an a priori unknown number of 

^'Bar' and 'bell' refer to image properties; we use 'bridge' and 'circle' to refer to the corresponding pieces of a dumbbell- 
shaped region. 
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Figure 10: Results on circle separation comparing the HOAC 'gas of circles' model to the classical 
active contour model. Top: the original images. The intensity of the bar takes values equally spaced 
between 48 and 128 from left to right; the background is 255; the bells are 0. In the middle: the best 
results obtained using the classical active contour model (8,1,1). Either the circles are not separated 
or the region vanishes. Bottom: the results using the HOAC 'gas of circles' model (2, 1,5, 4.0, 8, 8). 
All the circles are segmented correctly. 

circles of a certain radius, with short-range interactions amongst them. The requirement that circles 
be stable, i.e. local minima of the energy, fixes one of the prior parameters and constrains another. 

The 'gas of circles' model has many potential uses in computer vision and image processing. 
Combined with an appropriate likelihood, we have applied it to the extraction of tree crowns from 
aerial images. It performs better than simpler techniques, such as maximum likelihood and standard 
active contours. In particular, it is better able to separate trees that appear joined in the data than a 
classical active contour model. 

The model is not without its issues, however. The two most significant are related to the principal 
error modes found in the noise experiments of section 3.4: circles are found where the data does not 
ostensibly support them (false positives, a.k.a. 'phantom' circles), and two circles may be joined 
into a dumbbell shape and never separated. We discuss these in turn. 

The first issue is that of 'phantom' circles. Circles of radius fo are local minima of the prior 
energy. It is the effect of the data that converts such configurations into global minima. Were we 
able to find the global minimum of the energy, this would be fine. However, the fact that gradient 
descent finds only a local minimum can create problems in areas where the data does not support the 
existence of circles. This is because a circle, once formed during gradient descent, cannot disappear 
unless there is an image force acting on it. We thus find that circles can appear and remain even 
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though there is no data to support them. Adding a large level of noise exacerbates this problem, 
because random fluctuations may encourage the appearance of circles as intermediate states during 
gradient descent. 

The second issue is that of joined circles, discussed in section 3.5. Although the current HOAC 
model is better able to separate circles than a classical active contour, it still fails to do so in a number 
of cases, leaving a bridge between the circles. The issue here is a delicate balance between the 
parameters, which must be adjusted so that the sides of the bridge attract one another, thus breaking 
the bridge, and so that nearby circles repel one another at close range, so that the bridge does not 
re-form. Again, this is at least in part an algorithmic issue. Even if the two separated circles have 
a lower energy than the joined circles, separation may never be achieved due to a local minimum 
caused by the bridge. Again, high levels of noise encourage this behaviour by producing by chance 
image configurations that weakly support the existence of a bridge. 

We are currently working on solving both these problems through a more detailed theoretical 
analysis of the energy, and in particular the dependence of local minima on the parameters. 



A Details of stability computations 

In this appendix, starting from the equation for the circle and the expression for the radial perturba- 
tion in terms of its Fourier coefficients, 

-fit) = 70 W +(^7^ = {r{t),0{t)) = {ro{t)^Sr{t),Oo{t)) (A.la) 



where 



7o(t) = (ro(t),^oW) = (ro,t) (A.lb) 

and 

6r{t) = ^ake'^'^\ (A.lc) 

k 

with k G {m/ro : m G Z}, we give most of the steps involved in reaching the expression, equa- 
tion (2.6), for the expansion to second order of Eg around a circle. 
The derivative of 7 is given by 

d{t) = 1 (A.2a) 

r{t) = 5r{t) = ^akiroke'"^'^' . (A.2b) 

k 

The tangent vector field is given by 

r{t) = f(t)dr + e(t)de . (A.3) 
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We need the magnitude of this vector to second order. The metric in polar coordinates is given by 
ds'^ = dr'^ + r'^dO'^, so we have that |r(t)p = r(t)^ + r(t)^ by equation (A.2a). Substituting from 
equations (A.l) and (A.2b) gives 

|T(t)|2 = r2 + 2ro J2 + E «fe«fc'e'''°^'+''^*(l " r^o^k') . (A.4) 

k k^k' 

Taking the square root, expanding it as Vl + ^ + and keeping terms to second order 

in the ak then gives 

\T{t)\ = ro I 1 + ^ ^e*'-"'^* - \Y,auaykk'e'^'^^^+'''^' \ . (A.5) 

[ k k,k' J 

A.l Length 

Using equation (A.5), the boundary length is then given to second order by 

L(7)= £ dt\Tit)\=2nrol^l + ^ + l^e\ak\^^ , 
where we have used the fact that 

f dt e''""*'* = 2TrS{k) , (A.6) 

J —TT 

and that a-k = a^, where * indicates complex conjugation, because Sr is real. 

A.l Area 

We can write the interior area of the region as 



/TV Pr[U) PTT 

dO / dr'r' = / dO-r^{0) . 
-TT Jo J-TT 2 



Thus, using equations (A.l), and again using equation (A.6) to integrate Fourier basis elements, we 
have that 

A(7) = nrl + 27rroao + tt ^ | a/e P . (A.7) 



A.3 Quadratic energy 

To compute the expansion of the quadratic term in equation (1.3) for Eg, we need the expansions of 
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A.3.1 Inner product of tangent vectors 

The tangent vector is given by equation (A. 3), but we must take care as r(t) and T{t') live in dif- 
ferent tangent spaces, at 7(t) and 7(t') respectively. Since parallel transport does not preserve the 
coordinate basis vectors dr and de, it will change the components of r(t'), say, when we parallel 
transport it to the tangent space at 7(t). It is easiest to convert the tangent vectors to the Euclidean 
coordinate basis, 

dr = cos{0)dx + ^\ii{0)dy 

de = —r ^m.{6)dx + rcos{6)dy , 

as these basis vectors are preserved by parallel transport. Doing so, and then taking the inner product 
gives 

T • t' = cosiO' — 0) [rg + rodr + rodr' + 5r5r' + 5r5r ] 

+ shi(0' — 0)[rQ8r — tqSt + SrSr — SrSr'] . 

where unprimed quantities are evaluated at t and primed quantities at t\ Note that when t = t\ the 
expression reduces to equation (A.4). 

A.3.2 Distance between two points 

The squared distance between ^{t') and j{t) is given by 

\j{t')-j{t)\' = {x{t')-x{t)f^{y{t')-y{t^^^ 

= [(ro + 5r') cos((90 - (ro + Sr) cos{0)]'^ + [(ro + Sr') sm{0') - (ro + Sr) sm{0)]'^ , 

which after expansion gives 

hit') - ,it)\^ = 2rlil - cos(At))|l + -iSr + Sr') + ^.^(1 - cos(U) | ' 

where At = 0' — = t' — t. Expanding \/l -\- x + |x^to second order and collecting 

terms, we then find 

R{t,t') = |7(tO-7WI =2ro|sin(AV2)| + |sin(At/2)|((5r + ^rO + 4^(^^-^^0' , (A.9) 

4ro 

where A{z) = ^ ^-^^^^ 
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A.3.3 Interaction function 

Expanding ^{z) in a Taylor series to second order, and then substituting R{t^t') for z using the 
approximation in equation (A. 9), and keeping only terms up to second order in (^7 then gives 



At 



sm ■ 



where Xo = 2ro|sin(At/2)|. 
A.3.4 Combining terms 

Now let t') = r{t) ■ r{t')^{R{t^ t')). Combining the expressions already derived, we have 

G{t,t') = 

rl cos(At)$(Xo) 



Foo , even 



{5r + 5r') ro cos(At) |<l>(Xo) + ro 



At 



sm ■ 



Fio, even 



{Sr - Sr)rosm{At)^{Xo) 



Fii, odd 



+ {5r^ + Sr'^) ro cos(At) j ^A{At)^'{Xo) + ^ro sin^ (^) $"(Xo) 



At 



sm ■ 



F20 , even 



{SrSr') cos(At) <^ ^{Xq) + 2ro 



At 



sm - 



$'(Xo) - iroA(At)$'(Xo) + sin^ (^)$"(Xo)| 



+ {5r'5r — 6r6r) ro 



At 



sm - 



sin(At)$'(Xo) 



F22, odd 



((5r(5r - Sr'Sr) sin(At) < ^{Xq) + ro 



At 



sm ■ 



F2 1 , even 



F23 , odd 



(SrSr ) cos(At)$(Xo) . 



F24, even 
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where we have introduced the notation Fqo . . • i^24 for the functions appearing in the terms of G, 
and 'odd' and 'even' refer to parity under exchange of t and t' . Note that the F are functionals of <l>, 
and functions of ro and t' — t (but not t and t' separately). Note also that each line, and hence G, is 
symmetric in t and t' . 

The integral in the quadratic energy term is now given by //^^ dt dt' G{t^t'). We can now 
substitute the expressions for Sr and Sr in terms of their Fourier coefficients, Sr{t) = a/ee*^°^^ 
and Sr{t) = akiroke'^^^^^ . Due to the dependence of the F on t — t' only, the resulting integrals 
can be reduced, via a change of variables p = t' — t, to integrals over p. We note that in the terms 
involving Fio, Fii, F20, F22, and F23, the presence of the symmetric or antisymmetric factors in 6r 
and Sr' simply leads to a doubling of the value of the integral for one of the terms in these factors, 
due to the corresponding symmetry or antisymmetry of the F functions. For example. 



We therefore only need to evaluate one of these integrals for the relevant terms. Below we list the 
calculations for all the F integrals for completeness: 





which survives because Fqo is symmetric; 




k 




which survives because Fiq is symmetric; 
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// dtdt' 5r(t)Fii{t' -t)= // dtdt' "^akiroke'"^'^' Fii{t' -t) 

J J —TT J J —TT ^ 

= J2 (^kirok If dp dt' e^'-o'=(-P+*') Fn (p) 

/TT />7r 
dt' e'"'"'''' / dp e-'''°^P Fii{p) 

= Y,'^kirok2TTS{k) j dpe-'''°''P Fn{p) 
k 

= 0; 

j r dt dt' 5r^{t) F2o{t' -t)= If dt dt' ^ ^ afcafe'e^'"''('=+'=')* F2o(t' - t) 

k k' "^"^-^ 

k k' "^-^ 

= ^a/ea_/e27r / dp F2o{p) 

U J —TT 



k 



which survives because F20 is symmetric; 
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dtdt' 5r{t)5r{t')F2i{t' -t)= // dt dt' ^^auawe'^^'^^'^^'^'^ F2i{t' - t) 

= ^^akak' I r dpdt' e^'^oki-p+t')^irok't' ^^^^^^ 
k k' -^-^-^ 

dt' / dp €-'■•'0''^ F2i(p) 

k k' •'-^ 

u J —TV 



27r^|afc|2 f , 

U J —TT 



dpe-''-°''P F2i{p) ; 



// dt dt' Sr{t)Sr{t) F22{t' -t) = ff dtdt' Y,Y1 afeafe'^^ofce*'"«('=+'=')* F22{t' - t) 

fc fe' "'"'-'^ 

= E E akakdrok2TT5{k + k') f dp e--''(fe+fc')P _p22(p) 

fc fe' •'-'^ 
= 0, 

because with /c + /c' = from the delta function, the integral becomes one over F22 only, which 
vanishes due to the antisymmetry of F22 ; 
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/ / dt dt' Sr{t)Sr{t') F23{t' - t) = dt dt' CLkak'ir^k' e'^-'^^'^^'^'^ F23{t' - t) 

k k' J J— IT 

/TT r-TT 
dt' e''-o('=+'=')*' / dp e-'^'o'^P F23{p) 

= ^ ^afeafc.irofc'27r5(fc + k') ( dp e-"'°'=f 
k k' •'-^ 

= -27r^|afc|2irofc T dp e-'^--^^ F23{p) ; 
■^-'^ 

//" dt' Sr{t)5r{t') F2i{t' -*)=[[ dtdt'Y^Y^ akak'frlkk'e'''°^''^+'''^"> F2i{t' - t) 

k k' JJ-TT 

= -YYakak'rlkk' f dt' e'^'^^^^'^'' f dp e'^^'^P F24{p) 
k k' -^-^ -^-^ 

= -YYl Gikak'rlkk'2Ti8{k + k') f dp g-^^^^^ F24{p) 
k k' "^-^ 

= 2^^|a,|Vo'fe' r dpe-'^^^PF24{p) • 

Using these results then gives equation (2.5), which in combination with equations (2.3) and (2.4), 
gives equation (2.6). 
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